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NONLINEAR FOURIER ALGORITHM APELIED TO SOLVING EQUATIONS 
OF GRAVITATIONAL GAS DYNAMICS 

B.I. Kolosov 

Academy of Sciences USSR, Institute 
Space Research, Moscow 

Introduction 

The development of efficient numerical methods of solution 
of the equations of gravitational gas dynamics has- becomecparticu- 
larly urgent, in connection with the fact that the correct solu- 
tion of this_ problem models the behavior of a gas flow in the 
vicinity of 'a "black hole." In turn, this can be of direct 
interest in observational astrophysics. The strong nonlinearity 
of the problem complicates the possibility of detail'ed numerical 
investigation of the two dimensional flow pattern, and the neces- 
sity of introduction of artificial viscosity and the subsequent 
incorrect approximation of the boundary conditions in a certain 
finite vicinity of a gravitating center in a rough difference net- 
work leads to a qualitative difference in the results of the works 
ofwarious authors [1-3] . An attempt to overcome these difficul- 
ties was made in work [4] , where some self modeling solutions of 
the gas dynamics equations. were investigated, which could repre- 
sent the as 3 ?mptotics of flow in the vicinity of a gravitating 

center, on the assiamption of uniformity of the flow at infinity. 

In the present paper, the formulation of. a nonlinear algorithm vflth 
the use of a' Fourier series is given, which permits solution of 
the two dimensional stationary problem to be reduced to a system 
of common differential equations, with subsequent solution of the 
latter by computer. Within the framework of the approach developed 
below, a proof of ft.e existence of conical shock waves, with the 
cone vertex in the gravitating center, is successfully constructed' 
for the adiabatic index l<y<5/3. 

■Xd 

Numbers in the margin indicate pagination in the foreign text. 
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1. Formulation of the Problem 


Let, in region L.(0<r<«, e^t-Tr, tt] ) , axisyrametric , stationary 
gas flow with a gravitating center of mass M be described in a 
spherical coordinate system by the system of equations [4] 

-»'r Y 
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( 1 ) 


where the following notations are used: u^(r, 0), Ug(r, 0) are 

the radial and tangential components of the velocity; p(r, 9) is 
the gas density; P(r, 0 ) is the pressure; S'(r, e) is the entropy 
function; G is tie gravitational constant; r, e is the radius and 
angle of inclination of the radius vector to the axis of symmetry, 
respectively. It is assumed that functions p , S and P are con^,; 
nected by the equation of state 


/4 


P : 
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( 2 ) 


where y is the adiabatic index. 

It is required to find the solution of (l)-(2), with supple- 
mentary conditions at infinity 

, U© fe)=iio(e)st«e > 9 "- Pole), ^ ( 3 ) 

and the axial symmetry condition with 0=0 and 0=tt 


T> - C> 


'*'*+^ 1 . 
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(4) 


Further, instead of pressure P(r, 0), it is advisible to consider 
the enthalpy, after determining it by theerelationship 








(5) 


Now, in accordance with (2) and (5) , system (1) takes the form 


2 



( 6 ) 


E (r.9) - ^ 4 ^ ^ 

U» (iv . 

V ^ + vi v,1g Oi'o^e>)= o ^ 

»N^ 

+'h«3b? = o . 



where it is asstuned that 


tCr,9)^ VMM- 

r^ 

; Cv.e) - tiy -V Uq , ^ ~ 


QfM 

"T^ 

9CV,I3^ . » 


In concluding the formulation of time problem, we make the 
assumption that the solution of system (6) is smooth, namely, we 
will assume that the shock waves break down the entire region of 
determination L(r^[0, “],eg[-i:, it]) into £(£>1) subregions 1 ,^, 
within which t^ie functions u(r, e) , Ug(r, e), 6)? s(r, e) 

and H(r, e)‘ are infinitely differentiable with respect to angle 6, 



condition is valid 

C9Vn^-v “ > 

cp+9uiv . 



2. Formulation of Fourier Algorithm 

System of gravitational gas dynamics equations (1) and system 
(6) equivalent to it are systems of nonlinear equations, the numer- 
ical solution of which presents great difficulties. As is known 
[5] , to construct efficient algorithms for solution of nonlinear 
equations, the fortunate choice of the form of presentation of the 
desired function is of decisive importance. Therefore, we dwell in 
greater detail on the procedure of representation of the solution 
and its jiustif ication. For this purpose, we define velocity c> 
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components n^-Cr, 6 ) and' ..UgCr^S) by the expres-sion 

uv ct.e ) - - u ) co^ ct> 1 j 

UeC^iO)- ct>Cv;B) j 

where u(r, 0) is the velocity modulus, 0 ) is the angle of in- 

clination of u^(r, e) to vector u(r, e). 

'•'Wet reestablish the form of function 
§(r, 9). Available information in the 
formulation of the problem (l)-( 4 ) per- 
mits this to be done unambiguousjly . It 
fact, we analyze the geometric image of 
vector u. It follows from Fig. 1 that 

«i 5 tv,e) 5 : 0+ Cr,e). ( 9 ) 

Further, from boundary conditions ( 3 ) 
and S 37 mmetry conditions ( 4 ), the follow- 
ing properties of the functions flow 

- Eig. 1- 

I. ^ , 

2 » «P,Csr,0l= <D4t*-,0+5iTric>,j<=o,'i,.. , 

By virtue of the definition, angular function every- 

where unequivocal and finite (I^ij^Cr, e)|<A), with limited varia- 
tion of the function, piecewise differentiable in region L, with 
the exception of shock wave fronts «£-g(r, 63), on the boundaries of 
which, there are finite, unilateral derivatives 

"fc-^-i-0 -1-1’+ o — "I f 

In this manner, function ^]_(r, 0) satisfies the Lipshits continuity 
conditions [ 9 , 10 ] _ ^ — - 

d>Ur,0VOnCt.e>\^cU-e\ , oo,o^ci 4 -i, 
by virtue of which, there is the following proposition. 
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Theorem 1 


Function $][(r, 6), ever 3 rwhere in region L, excepting points 
0 = 0 i 2 , can be represented by the uniformly converging Fourier 
series 

Z MicClr'I'blYlkl?, (H) 

^ ^ [ 

and this representation is unique. ' 

In fact, conditions (10) satisfy only functions of the type 

and, since there is a countable set of such functions, (11) fol-.' 
lows from it, which proves the uniqueness of the representation. 

In turn, the Lipshits continuity conditions ensure uniform con- 
vergence of series (11) for e?^0g. We also propose that termwise 
differentiability of the series is possible. 

Further, by virtue of the uniform convergence, limited to 
the finite segment 

d>Cr,e)--e^ 2 I (12) 

for expression (9), we obtain the desired approximate Fourier /7 

representation 

■ ' cl 

UrC'f‘t0)=- UUr,0)CoS (0+ s I 

(.94 ^ J (13) 


3. Fourier Approximation of Two Dimensional Problem. (l)-(6) 
of System of Common Differential Equations 

We note that numerical solution of system of stationary gas 
dynamics equations (6) with a gravitating center assumes a transi- 
tion through the nonstationary problem of the method of establish- 
ment and difference approximation of the initial equations. By 
virtue of the fact that the equations to be solved have first order 
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differential operators, the difference equations which approximate 
them, because of limitations on stability, have the same first 
order of approximation* in difference; network . 

Fourier representation (11) -(13), established above, permits 
the, development of a fundamentally different approach to solution 
of the problem, which consists of reduction of the system of two 
dimensional equations in partial derivatives to a system of common 
differential equations in normal form which approximate them, solveji 
relative to the derivatives. Let the difference network with 
respect to angle 6 have Nq nodes. It is known [7] that any peri- 
odic function in such a difference network can be approximated by 
a discrete Fourier series, with number of terms N<N„ . Thus, the 
Fourier representation is the analog of'*a difference description of 
the problem with respect to angular variable 0 ."'Moreover , the 
uniform convergence of series (11) allows it to be hoped that the 
initial problem can be solved satisfactorily, with the minimum 
possible number of terms of representations (12)~(13). 

We also note a characteristic feature of equations (6) , which 
is that they include functions of different parity with respect to 
the angular variable, and that this circumstance leads to diffi- 
culties in obtaining a consistent order of approximation, with 
respect to^angle, of the approximating equations. In fact, it is 
easy to show that the application of an expansion in a Taylor 
series to system (6) gives fise to a meshing chain of equations. 

The mathematical technique developed below, within the framework 
of the Fourier representation under consideration, permits these 
difficulties to avoided and construction of a sufficiently general 
theory of reduction of gravitational gas dynamics equations (6) 
to a closed system of common differential equations, which can be 
solved with respect to the derivatives of the initial functions. 

We now proceed to description of the algorithm. For this 
purpose, We first substitute expression (8) in system (6), and we 
carry out identical transformations, with the use of the notations 
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( 14 ) 


I U^Cr.61 , 

I RlT,9)= ,, 

9 

where, following the work of Hunt [2] , the Bondi radius 0Rg=(^)^“X 
is adopted as the unit of distance, the asymptotic velocity of 
sound at infinity Cc is adopted as the unit of velocity and, in' 
this manner, we will assume further that system (6) is made dimen- 
sionless. 

Then, with (14) taken into account, from (6) there flows 



|t O 

I .cosd>Rir.0)= K.or.e), , 

* *% 

j- ^ 7 


Equation (19), in turn, is reduced to the form 

1 % *-v.K4H«r^ ril eg,' «r4; 

I - *- '«r ¥L'5f>i) ^ — 

- tos4)Hrr> 


(15) 

(16) 

(17) 

(18) 
(19) 


( 20 ) 


We introduce the notation 




T*©'* ’ 





'I 



( 21 ) 


Then, by solving equations (17) -(20) with respect to deriva- 
tive ■&$/3x, we ‘obtain 


A r = o , I 

j A (r.e) = A ir,e)' & (r,e ) , A tr.e) = ~ ,• | 

c (t.& j - ' H 0+ If a>) ■ . ' 


( 22 ) 



( 23 ) 




1 

_ Q-»i 


Here, we also write equations (15)) and (16) , which, with /_9 

notation (21) taken into account, take the form 




(24) 


As the following analysis shows, function s(r, 0) 


should be 


determined by tie relationship 



■XCte) 

E Ct :® } ’ 
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(25) 


where X(r, e) satisfies equation (24). 


Hunther, the algorithm consists of the following. We substi- 
tute representation (12) in equations (22)-(24). Now, with the 
use of symmetry conditions (4) , and the smoothness and infinite 
differentiability of the solution with respect to the angle every- 
where outside the boundaries of the shock wave front, we differ- 
entiate equations (22) -(24) in sequence, until, as 0^0 or 9->-ir, a 
closed system of common differential equations is obtained from 
variable r6(0, ») , for determination of coefficients m 2 h> 

^2u» 1, . . . , N, a=l, 2, N. For this purpose, the 

rule of Liebnitz [8] of differentiation of products of functions 
will be useful subsequently. It consists of the following 


where 




V --0 


_r'^\ /,< « 

- >-14) ) 


(26) 


Since the products of functions of different parities with 
respect to angle are included in equations (22) -(25), it is advis- 
ible to consider modifications of relationship (26) separately 
below, for the cases of even and odd derivatives. 
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1. Let u=2n+l, n=0, 1, and slm-ultaneously , Zj^(0)=OVk-2ji , 

&=0, 1, .... Then, 




\&rO- 


(27) 


2. Now, let u=X.n, n=0, 1, and )=01?lc=2il+l, £=0, 1, /lO 

. .. From this, we obtain 


^ in 
^ \fL V 9 =c 


(28) 


We proceed to calculation of the derivatives of function 
#(r, 6), for time being, bounded only by syiraietry conditions (4) 
at 0=0. By virtue of definition (12), we find 

<t>2wCo) =o, n = o,A,.. , 

+ C-^*' S. 7»^Crj 

Subsequently, the sequences of higher derivatives of functions 

O 

4=tg$ , F=l+tg $ , T=tg6 also are required, the values of which are 
determined by the following expressions 



I ^9-n o N, 

- (Pa , '^■5. «?5 =*Ps+ao£l>i<Pj+''& 

^6 - ^0<P,4>6+ e-4 <PA(i>bF4-t-A-^q)^ -V B ^ 

i *fi ~ ^ > 6 -j ~ H*l 6 j 

, ^ . ~i. — . . 

Now, since, with even derivatives of equation (22) with respect 
to angle at- 9=0, the identity satisfies 0, we examine the odd 
derivatives which, by virtue of (27) , take the recurrent form 


7 ! 


where 


( 31 ) 
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r- j ;jni X«) 

c,^ _ ^ I 


'ifc 


X - ~7 tJ‘ u r 




2-0l-8) 

101-0 - f,„ I 


>-t'J 


From equation (31) , 


/II 


follows, from which, in accordance with (22), we obtain 

ciWrH „ ^ 

( tls 


(32) 


(33) 


Ot^'j V, • 

L ,C.t.,„*8,„-'.f,r,„--i^ T»<^kVviVc. , 

(9 _ 3L ^ • f ' 

1*1 (jivuA 

Thus, it was shown that system of common differential equations 
(31) has a triangular matrix, which permits its trivial solution 
with respect to the derivatives of variable r of desired functions 

N . Similarly, a system is found for determina- 


m 


2n 


, n=0. 1, 


tion of moments E 2 


n 


ix 


“ fto ”=0,A,... tv/. 


(34) 


Finally, we proceed to determination of the coefficients of Fourier 

representation (12) in equation (23). Here, in accordance with 

(29) , we find ^ 

", ^ ^ A]h i«-i< 


K=i 


^ ^ 'll ^ 'UiaI 


• (35) 


^10 - C^)u ■ 




} ■ 


(36) 


Further, we write (35) in vector matrix form. For this purpose, 
we introduce the notations 
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/12 


^=C'®S®-V, 

r= M 6 . N / 

from which, finally, we obtain 



(37) 


where inverse matrix P~^ exists ever 3 rwhere, by virtue of the good 
dependence of P, and it is calculated once for the given order of 
approximation N of the problem. 


Thus, the' solution of two dimensional gravitational gas 
d 3 mamics equation (6) , in the case of random flow in regions of 
smoothness of the solution with respect to angle 6 , is reduced to 
system (3N4-1) of common differential equations in normal form, 
solved with respect to the derivatives of ir of functions m 2 ^, ^ 2 ^, 
u==0, 1, by means of which, at each point (r , 0)6i£'^t(r, 

e5, the two dimensional flow pattern can be regenerated 


UrCne)is- w <oiCB+ ^ i.mV'e ; 

. ^ ^ 

qe (yi® ) - vh ^ 0» e; ( e + .5: vb ) 


(38) 


To complete formulation of the algorithm for system (33) -(37), 
the initial values of the functions determined, which are produced, 
either by boundary conditions (3) atiiinf inity , or by the asymptotic 
conditions in the vicinity of the center (r^O) , must be added. 

¥e note that the problem under consideration is solved most simply, 
in the case of potential flow with uniform boundary conditions, 
when R(r, 0)=O and E(r, e)=E„,(0). Below, as an example which 
illustrates the approach stated, a system of unidimensional dif- 
ferential equations is presented, to which equations (6) are re- 
duced in the potential case, with N=l, /iful 





4\Sc- C< 


' / a -rti . 


(39) 
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“• tvvi4-c\)^ ^ 


Wlo- C< 


% 

C<i=0s^j(.-y ))'"oC«>)=M^^4MaCt^3-0^^.(;p,^„j5 

^ '1 * ^ 

- kftA^p/VT C\0&9ccm » ■ C-e. C««)= ''3 M ;i ,'j 


4. Existence of Shock Waves 


A characteristic feature of the gravitational gas dynamics 
equations under consideration is stationary shock waves, at the 
boundaries of which disruption of the continuity of the solution 
occurs. The nonlinear Fourier algorithm. (12)- (13), described above, 
permits analysis :and proof .of the existence of shock .waves by quite 
simple means. In fact, as has been pointed out, from the axial 
symmetry of the initial problem, the existence of two lines of 
symmetry flows, with 6=0 and 0=Tr. By using the latter situation, 
we now obtain from system (6) a system of differential equations 
approximating it, with the use of the e=ir line as the axis of 
symmetry. As a result, the same system of equations (33) -(37) 


is obtained, with the exception of 
in it , determined by relationships 


functions ^ 2 n+l~^ 2 n+l included 
(29) , wMoh should be substituted 


here by the expression 


i4 — ^ ^ a. IVH4 — 


(40) 


Now, let there be a solution of system (33) “(37) and continuity in 
region L, both With «&~(r, e, u~), and with $”^(r, 9, y”^) . But 
then, by virtue of the equality of the lines of symmetry at 0=0 
and that 6=tt, an angle 0g should exist, which divides the region 
of the determination of L into two subregions of smoothness L (r, 
0<0g) and L"*"(r, esTr-Og), at the boundaries of which joining of the 
solutions occurs from angular functions #^(r, 9, y”) and $'*'(r,e, 
y"*”) , which are not certainly equal to each other. In turn, angular 
functions andl#"^ are directly connected with the characteristics 
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of equations (6). Actually,' by definition, the characteristics 
(or flow lines) are expressed by the relationships. 

■ - — ^ ar _ Ao 

— 5 

i 

from which we obtain /14 

7^ ; ■ (41) 

Thus, each subregion L" and L"*" should be described by its family 
of characteristics r (0) and r (e), for which, in the limit of 
point r(9)=0, there is a common point of both families. Since, 
subsequently, analysis of the behavior of the solution in the 
vicinity of the center is of basic interest, we first note that 
there is asymptotic estimate for it [4] 

■= OCtf’j)' v-o, po. (42) 

On the other hand, by virtue of the fact that : 2 . }iuivfy*yiQ ^ 

is the change in angle of movement of the gas as a result of the 
gravitational force, which has to be finite everywhere, the exist- 
ence . of the finite limit 

©)| - \ ^ \ ^ oi <■<>>? , (43) 

follows from this. Finally, since, by definition, the shock wave 
front is connected with an abrupt change in the direction of move- 
ment of the gas, further, there is the following proposition about 
shock waves . 

Theorem 2 

Let, in region L(r^[0, ~] ,ee[— it, it]), with 1 iy^ 5/3, there be 
a solution of system of ax i symmetric, stationary, gravitational 
gas d 3 mamics equations (6) , which satisifies boundary conditions 
(3) at infinity and asymptotic behavior conditions (42) in the 
vicinity of the center 
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■wh'ep.e 


* y i- 

hi-s uC'ini', Ut:-vw ^cnOUa), U»-Vw'*‘V'V<t>Cir,9j f 

‘ 


Let further under the same conditions there be a solution of 
system of differential equations (33) -(37) whioh approximate them, 
for and on the lines of symmetry 0=0 and 9=^, respecr. 

tively, a 


... ^ 


} 

1 kf 1 ii 


u . 


Then, with w^(v) ^0Vk=2£+l , a curve 4s (r, 0)>O is found, which goes 
out of the gravitating center at angle of inclination 0g<ir, on 
which, with arbitrary N>1, angle there 

is the discontinuity |4>’^(r, this 

— * 4 “ 

case, in the regions L (e<9g) and L (9<0g), the solutions are 
described bycthe corresponding families of characteristics r (^q, 
9q, 0) and r"^(0, which has intersection points ry9g= 

r'l’Og on <£g(r, 9g) , which, by definition, means the existence of 

a stationary front of a conical shock wave, with the origin in the 
center r=0. 

We proceed to proof of the statements, on the basis of which 
we propose establishment of the impossibility of unbroken continu- 
ation of the asymfeotic solution with respect to angle 0 , from 
region L~ to region L"*". For this purpose, it is sufficient here 
to use the first equation of system (33) , which has the form 




( 44 ) 


where 


1 



Now, since, as r->0, as^mutnote (42) occurs, we obtain the 
estimates 


Y ko(v; W°,= 4 r-- 

^ W'o Wo 

= )=-sr. 


A 

'I ^ 


¥e 


introduce the notation “t" Q< Co)P’ 


Then, from equation (44), as r-^0, we find 

rv. ‘ 

!l!p 
-1 




(45) 


from which two characteristic solutions follow 


We 81 


2.- We ^ ^ 


(46) 


¥e initially investigate the first case. From the definition 

of 4>i , as r^O, there flows 

1 ’ , “ 

■|,A"wCo;w = -^^, 

V^r I 

Further, we come back to determination of” "angTes $(r, a) .d 

' s 6- + ^ JuJI fo) , 

1 kd» I ' ^ 



4>^C«) - 9- -Y (oJ '=a'HW e . 


/16 


By virtue of theorem 1, functions ^~(e) and $"^(6) are continuous 
everywhere in the region of determination 0 6[-ir, ir] . Into consid-’ 
eration of the function, we introduce 


and the notations 


, ' cb+(ej, e>0> 


I 

J 




, f- 2 c-u / k(o)k^ fr: 


(48) 


Then, it is evident that, by virtue of (37), there are the • 
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inequalities if only y^(0) ^ 0 for'yk=23L+l , £=0, 

1, We now show that function (48) is continuous ever3?where, 

with the exception of 6g6(-TT, "fr) • In fact, became of the contin- 
uity>'of $"(e) and ^"^(e) , for small angles, the expansion 

Olo)? e Ut i Mii(o)K)'^6(La/.) ec^i fj . 1 


is valid, from which there follows 

^^0 - ®sCC*^^S^)-:r'Sj. 4: 0 . 

In a similar manner, for the vicinity e=ir-0’, we obtain 

(?'4o) 4- 6 ; 

"f" 

Since functions are continuous by pairs everywhere it follows 
that function §(r, e), determined by expression (48), also is 

continuous everywhere, with the exception of angle 9g. Thus, it 
is has been shown that, in the vicinity of the gravitating center, 
the solution of system (33) -(37), which approximates the initial 
problem, cannot be unbrokenly continued by angle from one region 
of smoothness to the other and, consequently, for some angles, 
there is a discontinuous solution. In order to complete the 
analysis of this asymptote (mQ<2) , we-^return to equations (34) 
for the moments of energy. Since 


in the general case, 
sent at ion 


for E(r, 0), as r'>0, we obtain the repre- 


. 00)' 1 t , 




(49) 


We substitute (49) in system (34) . By virtue of the triangularity /17 
of the latter, to carry out its nontrivial solutions, fulfillment 

of the following conditions is necessary 

■ 

’ a ) 


p,- '‘’VC’I--.;, 


(50) 
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We now find index y. at which, together with •5’y=(y+l)/ (4(y-l)) , 
solvability condition (50) occurs. In fact, by definition, 

■ 6.y=- (5-3y) / (4(|-1) ) . On the other hand, it follows from (50) that 
0l=~l/2. From this,pwe find (5-3y) / (2 (y-1) )=1 or y=7/5=1.4. 

Thus, it has been established that characteristic solution 

■y 

(46) for mQf2 provides the required as;^ptotic behavior, only with 
index y=1.4. In this case, the flow in both regions of smooth- 
ness is nonpotential, and it cannot be continued unbrokenly from 
region L to L . Another confluence case is index y=5/3, at which 
6 u( 5/3)=1 and, consequently, This corresponds to potential 

fi'ow conditions in both regions, aiso separatediby a discontinuity 
boundary. 


We now proceed to the asymptotic solutions with m«=2, 6y#y+l 

“ V7-1) 

Here, for enthalpy H(r, 0i)), we obtain 

^ a ^ IT) " N P" 1,} 

: jpfi, = n:* 2, , (51) 


We note the asymptotic behavior of the pressure in the vicinity 
of the center was used as a boundary condition in works [1, 3] , in 
numerical solution of the two dimensional problem, where the pres- 
sure and, consequently, the enthalpy, in distinction from (51) , 
was erroneously assumed to equal zero (H(r<s, 6)=0. 


Further, in order to find the asymptotic behavior of function 
^k(O) ’ initials equation (24) for the energy. We note first 

that the solution of equation (34) and the approximate solution 
corresponding to it E(r, 8), approximated by finite Taylor series 
(49), exactly satisfies equation (24) at point 6=0 or 9'=Tr-0=O 
and approximately, at other values of angle 6 . In connection with 
this, we set the difference problem in the vicinity of the center, 
namely, we break down interval [0, 6g] into a series of nodes 0g/ 
(N-1) , where W is the order of approximation, Ogiir, and we demand 
that, at each value of angle 6 j=j 9 g/ (N-1) , approximate solution 
E(r, e) exactly satisfy equation (24). We then obtain 


( 
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or, with condition of nontrivial so-fvability (50) taken into 
account , 


where 


^ C«)‘ + 2 , (o)^v, V( Oj ^ ^ 

■ ] T ^ Ce'j + i t (0) ve] ) 


( 52 ) 


In this manner, in the vicinity of the center, the asymptotic 
values of coefficients u|^(0-),, k=l , 2, N are found from solu- 

tion of system of equations (52) . It is easy to see that u^(0)^ 
]ii^(0) . Actually, let N be sufficiently large. We assume that 
9j=9j and v^(0) =vi^(0) . But, since, by definition, <j> (® j » 

^'(6j^ Vtp and P2v^^u^^^2v^^k^ ’ follows from (52) that u"'"\0)^ 

Ic 

ij^CO) . Consequently, as above, the asymptotic solution cannot 
be continued without break from region L” to region L~*". Below, 
in the description of the shock wave front search algorithm, it 
will be determined that, in the vicinity of the center, there is 
a unique value of angle 6g, at which Hugoniot condition (7) is 
satisfied. 

We show now that each region of smoothness of solution L" and 
L"*” is described by its family of characteristics, the set of points 
of intersection of which lies on the associated ('conical) shock 
wave. For this purpose, we analyze characteristics (41) 

C.9-+ ) 

Since point r(o)=0 is a common point of family of characteristics 
(53) , here, for region L^, it is sufficient to determine 'the com- 
mon slope of characteristics r^, which go out from point r(0)=O. 

For this, in the vicinity r->-0, we consider angles 0=ir-e’, e'<l, 

~ =-^30k-oVi/i?^^v«(wH3‘j) = - -Iq (e>+ S, h? = 

QJl ^ k ' '' 

; = - ■*'3 C'l (“J )} ‘ 
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Since the flow is not potential in region L"*~, by virtiie of 
nontrivial solvability conditions (50), fn (iJi.(0)) = -l/2, it follows 
from this that dr>°^ i.e. , e.(r) decreases with increase in r. 

This corresponds to the slope of tie curves presented in region I 
in Fig . 2 . 



Fig. "2. 


We now make a qualitative 
analysis of the behavior of the first 
family of curves for large values of 
r(0). Because of the awkwardness of 
the expression, we limit outselves 
here to one term of the represen- 
tation of vi^(r), on the assumption, 
by virtue of conditions (3) , that 
it is small at large r and has an 
asymptotic form at infinity 


Then, 




^ 'I 




1 - 


from which, disregarding ^cos0 , we find 
where rQ(6Q)’is the sighting distance for angle 0g. 
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It follows from the resulting expression that r”(rQ, 0 q, 6) 
0:^0* Since, by virtue of estimate (45), y (r) has a finite limit 
in the vicinity r->0, here, we use asymptote y 2 ^(r)=ygl^, and we 
evaluate the behavior of the curves near the center. On the assimip- 
tion that 0=tt i , we obtain 




_ _ rsi'H 8^) 


Let yg>l. 
that 


It follows from this that a value of 0 ' 


, Ti-ev ^ ^0 




is found , such 
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which proves the increase of e'aand, together with it, 6(r)-_^+0' 
,as' r-J-O(y^). 

Thus, it has been established that there are two families of 
curves r“(rQ, 6 q, e) and r"'’(0, 6 q, 0), which move towards each 
other. Consequently, there must exist such angles 0g, that 

rj’Sg, which completes the proof of the theorem. 

Note 

Since the structure of equations (33) - (37) ipermits ^discontin- 
uous solution with respect to variable r, the existence of sta- 
tionary conical shock waves, with the apex in the center, does 
not exclude a frontal shock wave, and, thus, the pattern of gas 
flow from the gravitating center can be the superposition'! of two 
shock fronts. 


5. Conical Shock Wave Front Search Algorithm 


As was shown above, solution of the stationary gravitational 
gas dynamics equations in regions L, by means of the nonlinear 
Fourier algorithm, is reduced to the solution of approximating 
systems of common differential equations in subregions L“ and L , 
separated by shock wave boundary oi!g(r, 6g), the location of 
which was previously unknown. Let n and ? designate the unit nor- 
mal and tangent vectors at each point 0g). Then, con- 

dition (7) should be satisfied inig(r, 0g). Let, in a small 
•vicinity of the center, the shock wave front be described by the 
generatrix of a cone, which goes out from the center at angle of 
inclination 6g. Then, relationships (7), after simple fransfor- 
mations, are -written in the form 






d4)- 




(54) 
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( 54 ) 


where notations (14) and (25) are used 

E-(»ie>)= , e*V.0o-£e^(.j 

1 -hA) ^ 4 o S; ’ , 

X-A_, m-Vvi- , E=E-, H-Vl_ f 

± “ ’ ± 
d) 0^165) =: 9^4 v:0> 


(55) 
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Now, the algorithm of the solution of the complete problem 
briefly consists of the following. In regions L~ with boundary 
conditions (3) at infinity, the Cauchey problem (system (33) -(37)) 
is solved and, then, at each point (r, e), the approximate solu- 
tion is regenerated according to rejl:ationships (55) . We proceed 

-I- 

to regions L . Here, in the general case, the boundary conditions 

at infinity are unknown. Therefore, as the boundary conditions 

for the Cauchey problem in L^, asymptotic conditions (42) -(46) 

should be used, with the solution mQ=2 in the vicinity of the 

center. Let be bounded by angle 0g<^, which must be found. 

As the first step, we assume 6o=0 where is some randomly 

( 0 ) . ^ 

assigned angle ”^€(0, it). We designate 0g=7r-0g, and we break 

down interval [0, 9g] into (N-1) segments, so that 0j(6g)=j0'g/ 
(N-1) . Then, as was shown above, the asymptotic values of co- 
efficients ii^(0, 0g)aare found from solution of system (52) .an 

i V (56) 

, N-/1 

In turn, from relationships (54), for functions B (9g)=c0ss^(6g) , 
we obtain an expression which, with the arithmetic calculations 
omitted, is written in the form 

hr- a , ^ ^ , (57) 

i ( tsi> [ cos d> ’ 


y 
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Where the following designations are used 


U + W d) j ?■='?>- * ' 

<1 1»0 - Jua.aaC^,).^C^-|g;,g0^^oW ' 


J21 


(58) 


Since, in passage of the shock wave front, the inequality 
m_^<m must occur, p( 0 g)>cosf is equivalent. Therefore, further, 
we assume condition g^l is satisfied simultaneously 

here,. Now, since the values of sre known from solution 

of system (56), with given equality 

we find a new value of If it turns out that 6g h6g t we 

B {^1') 

assume 63=63 ^ and, again, we find the solution of system (56) and 

Cl) 

3 ( 6 g We will continue this procedure, until the desired ^alue 

( o\ 

0 g=lim 0 g^ ' is found. 

The next step is the search for asjzmptotic values of m^_(r, e) , 
E_|_(r, e) and \^(r, 0 ) in the vicinity of the center. For this 
purpose, it advisable to return to initial differential equations 
(15)~(20), which, with the asymptotic behavior of the solution 
taken account, take the form 

+ (60) 

Xe ’=■ 

j gu X0 m-C5-«)>i ^ ^ ' 

where, in reducing requations (18) ^( 20 ) to equations (60) - (61) , the 
properties of equations ( 59 ) were used, and e was a previously 
assigned small number/ _ 

It can be shown that, from sjrmmetry requirements ^-j =0, H^=0 

6= S 

follows . 
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Thus, the problem in the vicinity of the center was determined 
in the- general case, to within random Constants Eq and Xq, for the 
selection of which, the involvment of some physical considerations 
if required. By virtue of the continuity of the energy on the lines 
of S 5 nnmetry, we subsequently will assume that Eq=Eq, >>q=Xq. As /23 

analysis shows, system (59) -(61) is solved nxmerically more simply 
in order to obtain the values of the coefficients sought in expres- 
sions (55) fr.om the numerical results. Therefore, here, we break 
down interval [0g, tt] into N angular segments (Tr-6g)/N, so that 


0j=0g+jAe. Then, by using the scheme of second order of approxi- 
mation of equations (59) , we obtain the difference scheme 

■ ^ , (62) 

for the solution of which, in turn, we use the trial run method 
17] 


where 


Ejl - QtcE||t+/i -V It, , 


- 0 - 


V- 






ao=o,i„=E_0r,ei) 


(63) 


The solvability of (63) , as is known, is connected with stability 
of the trial runs, which always occurs with 0g©(Tr/2, it), since, in 
this case, (9j^)<0, k=l, 2, ..., N. Equation (59) is 

solved for X(r<e,0), by means of the same algorithm. After the 
numerical values of Ej^=E(r<e, e^) , Xj^=x(r<e, 6^^) are found, we 
approximate them by series (55) , and we then determinetthe values 

The asymptotic values of m_^(r, e) and H_j_(r, 6) are found in 

a similar manner where relationships (54) are used as the boundary 

conditions for equations ‘('60) -(61) at 0=0c; and, with e^-ir , from 

^ 0 

the asymptote we have m^(r<e, ■ir)=2, H_^(r<e, iir)=HQ=0. Here, 

"V 

in distinctionffrom equation (59) , together with the requirement 
for second order of approximation, an iteration scheme for solva- 
bility of the difference analog of equations (60) -(61) also should 
be used. 
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Thus, the solutions of equations (59) -(61)^ completely define 
the problem in region in the vicinity of the center and are the 

initial conditions for solution of system of equations (33) -(37). 

Since the shape of the surface of the shock wave front can change 
with distance from the center, it is advisible further to deter- 
mine the position .of the shock wave from the intersection of the 
r (0, 6q,0) and r< C^q, 0q, 6) curves which, in turn, permits 
determination of the flow in the total region L(0<r<“, 6£[-tt, it]). 

The algorithm reported above was programmed for a BESM-6 /24 

computer. The results of the calculations and discussion of them 
will be published separately. 

In conclusion, I thank G.S. Bisnovatyy-Kogan, who drew the 
attention of the author to the characteristic features of the 
problem, an well as to V.A. Ladygin, for discussion of some ques- 
tions . 
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